# -*- coding: utf-8 -*-
"""
Created on Mon Oct  8 14:32:52 2018

@author: luolei

CSTR反应器模拟
"""
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import pandas as pd
import numpy as np
import sys

sys.path.append('../..')

from lib.cstr_simulation import *

seg_len_range = [30, 50]  # **手动设置每次运行时间步数变动范围


def random_seg_len():
	"""生成每次平稳操作运行数据长度"""
	seg_len = np.random.randint(seg_len_range[0], seg_len_range[1])
	return seg_len


def random_op_params(op_params):
	"""随机改变一个操作参数"""
	chosen_loc = np.random.randint(0, operable_n)
	chosen_op_col = op_cols[chosen_loc]
	col_bounds = op_col_bounds[chosen_op_col]
	op_params[chosen_loc] = np.random.uniform(col_bounds[0], col_bounds[1])
	return op_params
	

if __name__ == '__main__':
	#%% 初始状态
	ca, T = ca_init, T_init     # 初始状态参数：反应浓度和温度
	op_params = [ca_0, T_0, q]  # 初始操作参数：进口浓度、冷却水温度、冷却水流量

	#%% 设备运行
	op_param_records, output_records = None, None
	while True:
		# 设定每次平稳操作参数运行时间长度
		seg_len = random_seg_len()
		t = np.arange(0, seg_len * dt, dt)

		# 积分求解输出变化
		output_seg = odeint(anisothermal_reaction, (ca, T), t, (op_params,))
		output_seg = output_seg[:seg_len, :]
		
		#%% 数据记录
		# 操作参数记录
		op_param_seg = np.array(op_params * seg_len).reshape(-1, operable_n)
		if op_param_records is None:
			op_param_records = op_param_seg
		else:
			op_param_records = np.vstack((op_param_records, op_param_seg))

		# 模型输出变量记录
		if output_records is None:
			output_records = output_seg
		else:
			output_records = np.vstack((output_records, output_seg))

		# %% 更新状态变量和操作变量
		# 更新变量初始值作为下一次模拟的起始状态
		[ca, T] = list(output_seg[-1, :])

		# 更新操作参数
		op_params = random_op_params(op_params)

		#%% 终止条件
		if output_records.shape[0] > steps:
			output_records = output_records[: steps, :]
			op_param_records = op_param_records[: steps, :]
			break

	#%% 绘制参数和变量图
	total_n = operable_n + target_n
	loc = 'upper right'
	plt.figure(figsize = [6, 1 * total_n])

	# 操作参数画图
	for i in range(operable_n):
		plt.subplot(total_n, 1, i + 1)
		plt.plot(op_param_records[:, i], 'r', linewidth = 1.0)
		plt.legend([op_cols[i]], loc = loc, fontsize = 8)
		plt.xticks(fontsize = 8)
		plt.yticks(fontsize = 8)

	# 状态变量画图
	for i in range(target_n):
		plt.subplot(total_n, 1, i + operable_n + 1)
		plt.plot(output_records[:, i], 'b', linewidth = 1.0)
		plt.legend([target_cols[i]], loc = loc, fontsize = 8)
		plt.xticks(fontsize = 8)
		plt.yticks(fontsize = 8)

		if i == target_n - 1:
			plt.xlabel('time', fontsize = 8)

	plt.tight_layout()
	plt.savefig('../../graph/process_vars.png', dpi = 450)

	# %% 保存数据
	data = np.hstack((op_param_records, output_records))
	data = pd.DataFrame(data, columns = op_cols + target_cols)

	data.to_csv('../../data/runtime/cstr_data.csv', index = False)
	
	#%% 归一化处理
	all_cols = set(target_cols + selected_cols)
	data_nmlzd = data.copy()
	for col in all_cols:
		min_value, max_value = var_bounds[col][0], var_bounds[col][1]
		print('min {}: {:.4f}, max {}: {:.4f}'.format(col, min_value, col, max_value))
		normalize = lambda x: (x - min_value) / (max_value - min_value)
		data_nmlzd.loc[:, col] = data_nmlzd.loc[:, col].apply(normalize)
	
	data_nmlzd.to_csv('../../data/runtime/data_nmlzd.csv', index = False)
	
		
	
	
	


